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Abstract 

We have carried out a thorough benchmark of the FDE-ET method for calculating hole 
transfer couplings. We have considered 10 exchange-correlation functionals, 3 non-additive 
kinetic energy functionals and 3 basis sets. Overall, we conclude that with a 7% mean relative 
unsigned error, the PBE functional coupled with the PW91k non-additive Kinetic energy 
functional and a TZP basis set constitutes the most stable, and accurate level of theory for 
hole-transfer coupling calculations. The FDE-ET method is found to be an excellent tool 
for computing diabatic couplings for hole transfer reactions. 
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1 Introduction 


The quantum mechanical study of realistically sized molecular systems has become a goal for 
quantum chemistry and material science. To that aim, multilevel and multiscale algorithms 
(such as QM/MM) generally approach the problem by representing the system as a set of 
subsystems whose interaction is accounted for approximately. Along these lines, the frozen 
density-embedding (FDE, hereafter) formalism developed by Wesolowski and Warshel^ (see 
Ref. El for a recent review), has become a popular avenue of research. FDE has been ap¬ 
plied to a vast array of chemical problems, for instance, solvent effects on different types 
of spectroscopy,®^ magnetic properties,®^ excited states,^^®^ charge transfer states.®®^ 
Computationally, FDE is available for molecular systems in ADF,®® Dalton,l2^J^ and Tur- 
bomole^^SESI packages, as well as for molecular periodic systems in CP 2 K®* 2 ZI and fully pe¬ 
riodic systems (although in different flavors) in CASTEP,®^ Quantum Espresso,®® and 
Abinit.®® 

The FDE method casts itself in the framework of subsystem density functional theory, by 
which the electron density of the total system is split into subsystem contributions and can be 
determined by solving coupled equations featuring an effective embedding potential. In this 
way, polarization given by the interaction of the subsystems is included. This subdivision 
of the total electron density into subsystem contributions has lead to the use of FDE as an 
effective charge and spin localization techniqueAlthough the reasons for the ability 
of FDE to yield charge localized states will be addressed in the following section, here we 
will take this for granted and discuss why such a property of this method is of interest. 

The quest for computing charge localized electronic states has a long history,®® es¬ 
pecially in recent years with the advent of the Generalized Mulliken-Hush method® and 
constrained DPT.®® Charge localized states are also known as diabatic states because it 
has been shown that they minimize the corresponding non-adiabatic coupling matrix®® - 
a defining property of diabatic states.®® Modeling of charge transfer (CT) reactions often 
involves the use of only two diabatic states: a state where the charge is on the donor (D) 
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also called initial state, and a state where the charge is on the acceptor (A), also known 
as hnal state. However, because the charge localized states are not necessarily eigenstates 
of the molecular Hamiltonian of the system, nor they are constructed to be orthogonal, we 
expect the Hamiltonian and overlap matrix to be non-diagonal. Namely, 


H 


/ \ 

Hdd Hoa 


Had Haa 


, s = 


1 Sda 
Sad 1 


( 1 ) 


Whether the CT rate is computed with Marcus theory® or it is extracted from a non- 
adiabatic dynamics,®® it is related to the following matrix element:®® 


Vda ^ (gcA - 2 ■ (2) 

which is known as transfer integral, coupling matrix element, charge transfer coupling, etc. 
For many systems of interest, Vda depends strongly on the molecular geometry,®® thus 
the dynamic charge transfer process is more straightforwardly modeled with real-time dy¬ 
namics methods, such as Tully’s surface hopping or Ehrenfest dynamics® as advocated by 
many.®® For such a model to be computationally efficient, the electronic couplings between 
the diabatic states have to be computed efficiently. Several methods have been proposed for 
this task, such as semiempirical methods,®® methods exploting the frozen orbital approxi¬ 
mation®® (i.e., the molecular orbitals of the diabatic states are approximated by the ones of 
isolated donor and acceptor fragments), or all-electron methods such as wavefunction meth¬ 
ods which are subsequently rotated to yield diabatic states,®®® or the ones that focus on 
constructing the diabatic states by imposing locality of the electronic structure.®®®® 
This explosion of methods for calculating the coupling matrix elements between diabatic 
states for electron transfer processes called for the setup of a benchmark set which can be used 
by all researchers developing novel algorithms. In addition, a question which is frequently 
posed is whether charge transfer couplings are sensitive to the particular method chosen 
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for their evaluation. These questions were posed to an audience at the 2011 conference on 
“Charge Transfer in Biosystems” organized by Profs. Rosa di Felice and Marcus Elstner. 
The seed planted in that conference gave rise to an interesting work by Kubas et al.^ 
where the authors compared benchmark values of coupling elements for 11 molecular dyads 
calculated with correlated wavefunction methods associated with the Generalized Mulliken- 
Hush diabatization with values computed with the constrained density functional (CDFT) 
method, fragment-orbital DFT (FODFT) and density functional tight-binding (FODFTB). 
The test set was named HABll, and it was found that the all-electron CDFT method (as 
implemented by Oberhofer et al^ in the CPMD software^SSEOl) reach a mean relative 
unsigned error (MRUE) of 5% if it is employed in conjunction with 50% of Hartree-Fock (HF) 
exchange in the PBE functional, and a deviation of about 39% if the pure PBE functional 
was used. The non-selfconsistent fragment orbital method yielded a deviation of 38% and 
the semiempirical FODFTB of 42%. 

HABll consists of eleven 7r-conjugated dimers, plus four additional aromatic rings. In 
Table [T] the structures of the monomers are shown. These organic molecules were chosen 
because they feature different vr bond arrangements and different kinds of heteroatoms. The 
monomers are: ethylene, acetylene, cyclopropene (having one high electronic density bond), 
the antiaromatic ring cyclobutadiene, O, N and S containing heterocycles, hve polycyclic 
aromatic hydrocarbons (benzene, naphthalene, anthracene, tetracene and pentacene), and 
one derivative of benzene: phenol. These organic compounds are well known to be part of ef- 
hcient semiconductor materials^^^®^ and some take part in CT processes in biomolecules.*^^*^ 
We obtained the Cartesian coordinates for every structure from the reference^ (for details 
about how the geometries were obtained we refer the reader to that source). 

The purpose of this work is to provide the community with information about the per¬ 
formance of the FDE method against the HABll benchmark set. As the FDE method is 
used to obtain the electronic structure of the diabats, a post-SCF calculation follows for 
the determination of the couplings. l^^*^^*^^ l The resulting composite method is termed FDE- 
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ET, hereafter.^ As we will see in the results section, the FDE-ET method compares quite 
well with the benchmark calculations, albeit experiencing outliers, in most part traceable 
back to convergence issues. To offer as complete of a picture as we can, the FDE coupling 
calculations are carried out with 10 different exchange correlation density functional(XC, 
hereafter), 3 non-additive Kinetic energy functionals (NAKE, hereafter) and 3 basis sets. 
This totals to a staggering 90 levels of theory tested in this work, leading to a total of 5400 
coupling calculations for the HABll set alone. In addition, and following Ref. [63, we have 
carried out calculations for dyads whose monomers were rotated with respect to each other, 
presenting an additional 3780 coupling calculations. 

This paper is organized as follows, in section]^ we show briefly the characteristics of FDE- 
ET. In section the computational details are described. Section collects the results of 
the comparisons against the HABll test set and for rotated ethylene and thiophene dimers. 
Finally, in section we outline the conclusions. 

2 Diabatic states from Frozen Density Embedding 

2.1 Background on FDE 

The FDE formalism prescribes the total electron density to be expressed as the sum of 
subsystem electron densities.l^^^^^ Namely, 

Ns 

Ptot(r) = ^p/(r). (3) 

1=1 

Where Ng is the number of subsystems. 

The electron density of each subsystem is obtained by solving a Kohn-Sham (KS) like 
equation augmented by an embedding potential that accounts for the interactions of other 
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Table 1: Dimers of the HABll test set. 


Dimer 

Ethylene 

Acetylene 

Cyclopropene 

Cyclbutadiene 

Cyclopentadiene 

Furane 

Pyrrole 

Thiophene 

Imidazole 

Phenol 

Benzene 

Naphthalene 

Anthracene 

Tetracene 

Pentacene 


SymboP 

(EE) 

(AC) 

(CP) 

(CB) 

(CD) 

(FF) 

(PY) 

(TH) 

(IM) 

(PH) 

(BB) 

(NN) 

(AA) 

(TT) 

(PP) 


Structure 
H2C = CH2 

HC=CH 



OH 

6 

o 



Reference method^ 
MRCI+Q 
MRCI+Q 

MRCI+Q 

MRCI+Q 

MRCI+Q 

MRCI+Q 

MRCI+Q 

NEVPT2 

NEVPT2 

NEVPT2 

NEVPT2 

SCS-CC2 

SCS-CC2 

SCS-CC2 

SCS-CC2 


a Abbreviations used in this work. 
b Ref. EZl 
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subsystems whose density is kept frozen in this step, such as 


-V 


V 


KS 


(r) + n, 


I ( 

emb \ 


= e(i)7(r)0p)7(r). 


( 4 ) 


Where are the molecular orbitals of subsystem I, and is the embedding 

potential acting on the same subsystem reading as follows: 
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I 

emb 



pj 
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+ 


5T,[p] 6T,[pj] 5E,,[p] 5 EM 

6p{r) 6pi{r) 6p{r) Spi{r) 


( 5 ) 


In the above, Tg, and are kinetic and exchange-correlation energy functionals, and 
the nuclear charge, respectively. A special comment for the kinetic energy is needed. In 
the KS method, Ts[p] should be calculated from the molecular orbitals of the entire system. 
However, these orbitals are not calculated in FDE and therefore are not accessible. Approx¬ 
imate kinetic energy functionals are employed instead, representing the NAKE term with 
a semilocal functional. This approximation is ultimately the biggest difference between an 
FDE and a full KS-DFT calculation of the supersystem.l^®^!] por example, when the sub¬ 
systems feature a large overlap between their electron densities, FDE in conjunction with 
GGA NAKE functionals becomes inaccurate when compared to regular KS-DFT.I ^^ * ^^ * ^^ I To 
achieve selfconsistency, the subsystem densities are determined in an iterative way called 
freeze-and-thaw.l^^*^ 


2.2 How does FDE generate diabats? 

Diabatic states can be generated with FDE by construction. In practical terms, the calcu¬ 
lation is performed on at least two interacting subsystems (donor and acceptor fragments) 
whose electron densities are determined via the freeze-and-thaw procedure employing ap¬ 
proximated NAKE functionals. A set of two simulations is set up: one featuring a hole 
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on the donor fragment (i.e., the KS-like equations are solved imposing the density of the 
corresponding subsystem to integrate to a number of electrons defecting by one compared 
to the neutral fragment), and another calculation in which the acceptor is now positively 
charged. We should note that it is also possible to increase the number of electrons by one - 
in this case, an excess electron is generated on the subsystem rather than a hole. The result 
of the calculations is that the hole (or electron) is completely localized onto the fragment it 
was placed on at the onset of the calculation. 

There are four reasons for the FDE calculations to yield charge localize states.^ First, 
the subsystem orbitals are not imposed to be orthogonal to orbitals of the other subsystems. 
This is important as it implies that not imposing orthogonality removes a bias towards de¬ 
localization, as noted by Dulak and Wesolowski.l^Sl However, this reason alone is not enough. 
A second reason is the fact that FDE calculations are carried out in the monomer basis set 
[i.e., using the FDE(m) method^. With no basis functions on the surrounding frozen sub¬ 
systems, a charge transfer between the subsystems becomes an unlikely event and the SCF is 
biased to converge to a charge localized solution. The third reason is similar to the previous 
one and invokes the fact that FDE calculations are always initiated with a subsystem local¬ 
ized guess density. The initial conditions also have a bias effect on the hnal SCF solution - a 
localized initial guess density will likely yield an SCF solution that is subsystem localized as 
well. The fourth reason is more subtle. It deals with the shape of the embedding potential 
in the region of the surrounding fragments. Electrons remain localized also because there 
are repulsive walls in the vicinity of the atomic shells of atoms belonging to the surround¬ 
ing subsystems. As noted by Jacob et a/.,®'the approximate kinetic energy functionals are 
unable to cancel out the attractive potential due to the nuclear charge in the vicinity of the 
nucleus. However, the shape of most semilocal kinetic energy potentials is such that in going 
towards the nucleus they start out too low compared to the exact potentials, then cross the 
exact potential and become larger in the region of an atomic shell. This is so up to when 
the shell has faded, then the potential becomes again too attractive (see for example Figures 


4 and 5 of Ref. IHZIand Figure for a simplified depiction of this effect). This behavior was 
also reported by Fux et when they calculated the approximate vs. accurate potentials 
for selected dyads. 

Until now, we have assumed that the the SCF procedure in FDE (i.e. when the so-called 
freeze-and-thaw cycles are used, see the computational detals) leads to a unique solution. 
This has been challenged recently in the limit of exact NAKE functionals.^ In this work, 
however, we only employ approximate functionals for which experience shows that a unique 
solution is always found. Numerical evidence of this can be found in Ref.®' 




emb 





Atom in a Frozen Subsystem 

Figure 1: Depiction of the shape of the embedding potential (in red) in the region of the 
atomic shells of surrounding subsystems (aka frozen subsystems). 


2.3 Coupling calculations with FDE: the FDE-ET method 

In case of non-orthogonal functions, the coupling (exact at the Hartree-Fock level but only 
approximate at the DFT level) can be expressed as:'^^®®' 

HnA = {'ipDmA) = SdaE • ( 6 ) 

where H is the molecular electronic Hamiltonian, 'ipD and 'ijjA are the two diabatic states (D 
for donor, A for acceptor) and p*'^^^(r) is the transition density which reads as: 

electrons 

= I(iPd\ ^ 6{rk - r)\'iljA) (7) 

k=l 
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Assuming the wavefunctions to be expressed in terms of single Slater determinants, the 
overlap element appearing in Eq. Q is found by computing the following determinant: 

Sda = det , (8) 

where is the transition overlap matrix in terms of the occupied orbitals 

Thus, the transition density is now written in the basis of all occupied orbitals 
which make up the diabatic states 'ipD and 'll! a- Namely, 

occ 

kl 

Finally, the Hamiltonian coupling is calculated by plugging Eq. @ and Eq. ([^ in Eq. and 
the resulting matrix elements in Eq. (|^ - that is, the coupling of two Lowdin orthogonalized 
V’d and In the FDE-ET method, the and 0^"^^ orbitals are borrowed from the 

FDE subsystem orbitals as prescribed by Refs.^^^® 

3 Computational Details 

All effective couplings were calculated with the Amsterdam Density Functional (ADF) pack- 
g^ggESl (2014 release). For this benchmark study, we selected ten XC functionals as follows: 
three GGAs (PBE,I2SI BLYP™] and PW9llS21), two MetaGGAs (M06-L™ and TPSSflS!), 
three Hybrid functionals (B3LYP,I^ BHandHH® and PBEQI^^ and two MetaHybrid func¬ 
tionals (M06-2X^I® and MOb-HF^I^. These functionals were employed in the FDE calcu¬ 
lations. However, the non-additive term for the non-additive exchange-correlation energy 
functional needed for the embedding potential and functional is computed always 

at the local or semilocal level®^ in order to avoid costly OEP type procedures. Specihcally, 
when B3LYP and BHandH were used, we chose BLYP for the PBEO was replaced 

by PBE and for both MetaGGA and MetaHybrids, PW91 was chosen. 
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Three NAKE functionals were used, the LDA Thomas-Fermi^ functional (TF, here¬ 
after), the GGA functional PW91k*^ and the gradient expansion approximation (GEA) 
functional P92.l^ 

Regarding basis sets, we use three Slater-type basis sets, TZP, TZ2P and QZ4P. TZP 
and TZ2P are double ( in the core and triple ( in the valence whereas QZ4P is triple ( in 
core and quadrupole ( at the valence. Additionally these basis sets are augmented with one, 
two and four polarization functions respectively.*^ 

The ADF default settings for the self-consistent field (SGF) cycles procedure were used. 
Also, the Becke*^ numerical integration grid and the ZlmFit*^ density htting options were 
set for both FDE and FDE-ET (through the ElectronTransfer keyword) calculations. 

The FDE-ET electronic couplings are obtained hrst by running an FDE calculation [i.e., 
solving for Eq. (|^], where the density is minimized by three freeze and thaw cycles, the hrst 
one of these cycles was performed using Thomas-Fermi NAKE and the next two were carried 
out with the corresponding NAKE. Secondly, a post-SGF evaluation of Eq. ([^ where each 
term is given by the 2x2 Hamiltonian and overlap matrices of Eq. ([^. In the FDE-ET 
post-SGF step, hybrids, metaGGAs, and metahybrid are not yet supported. Hence, the XG 
functional was changed in the evaluation of Eq. ([^ by a GGA functional. Specihcally: for 
B3LYP and BHandH the BLYP functional was used. The PBE functional was used when 
PBEO, MetaGGAs and MetaHybrids were employed. 

All dimer structures were taken from the HABll databases in Kubas et a/.*^In total, 15 
TT-systems perfectly stacked were analyzed, distance dependence calculations of the electronic 
coupling at 3.50, 4.00, 4.50 and 5.00 A were performed. In addition, we took ethylene (EE) 
dimer and thiophene (TH) dimer to compute the dependence of the coupling for different 
rotations with respect to the center of mass of each monomer in the same way as Kubas et 

a/.EZI 
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4 Results and Discussion 


The FDE-ET couplings are compared for each system against correlated ab-initio wavefunc- 
tion methods obtained by Kubas et al^ In that work, MRCI+Q, NEVPT2 and SCS-CC2 
were employed depending on the system size. MRCI-I-Q was used for ethylene (EE), acety¬ 
lene (AC), cyclopropene (CP), cyclobutadiene (CB), cyclopentadiene (CD) and furane (FF) 
dimers. NEVPT2 was used for pyrrole (PY), thiophene (TH), imidazole (IM), benzene (BB) 
[]and phenol (PH); hnally, SCS-CC2 was employed for the larger rings, such as naphthalene 
(NN), anthracene (AA), tetracene (TT), and pentacene (PP). For sake of completeness, in 
the supplementary information^^ we report all calculated couplings, all correlation plots 
and all the error analyses computed for each dimer. There, we report the mean unsigned 
error (MUE) and mean relative unsigned error (MRUE), the mean relative signed error 
(MRSE), and the maximum unsigned error (MAX) varying each of the following categories: 
XC functionals, the NAKE functionals, and the basis sets. To aid our explanation of the 
results, we have chosen to report in the hgures of the main text the variation of the three 
categories from a common starting point: the PBE/PW91k/TZP level of theory (e.g., XC 
functional/NAKE/basis set). 

4.1 Effect of the non-additive Kinetic energy functional 

Let us hrst discuss the behavior of the NAKE functional. In Figures[^(|^, the MUE (MRUE) 
for couplings obtained using the PBE functional and the TZP basis set are reported. All 
other data are reported in the supplementary materials.^^ A very good relation with the 
reference data is clear. The MRUEs are always below 20% and generally gets better and 
better as the distance separating the monomers increases. 

Overall, the NAKE functionals are consistent with each other. Each NAKE features 
couplings that correctly decay exponentially with the inter-monomer distance. For some 

Additional reference data for benzene using SCS-CC2 method is given in reference [STl In this work we 
compared with NEVPT2 level of theory. 
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Figure 2: MUE values as a function of different NAKEs used in this study. In this plot, the 
PBE functional and TZP basis set are employed. Full results are available in the supple¬ 
mentary information section. All bars are in meV. 
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3.50 4.00 4.50 SA 


3.50 4.00 4.50 5.00 


3.50 4.00 4.50 5.00 


Figure 3: MRUE in the performance of each NAKE functional. The PBE functional and 
TZP basis set are employed. Full results are available in the supplementary materials section. 
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functionals, we found a not so good description of the electronic coupling for a few dimers 
(e.g., AC, TT and BB systems in Figure S2 in the supplementary materials), we can see that 
again the NAKE functionals are consistent and do not change the picture in going from one 
NAKE to another. Nevertheless, in systems like TH or TT the Thomas-Fermi functional 
performs poorly when used in conjunctions with several XC functionals (Figure S2). We can 
attribute this behavior to the fact that the Thomas-Fermi potential compared to the GGA 
NAKE potentials is too soft and is not successful in localizing the hole, especially when the 
QZ4P basis set is employed. 

4.2 Effect of the basis set 

We now discuss the effect of the basis set in FDE-ET coupling calculations. Once again. Fig¬ 
ures |^([^ report the MUE (MRUE) when varying the basis set employing the PBE/PW91k 
functionals. The picture is not very different from the previous section in the sense that all 
MRUEs are at or below 20% with an inclination for being larger for aromatic dimers. Some 
improvements in the electronic couplings at long distances are noticed when the number of 
basis functions are increased. This can be explained by the fact that the more diffuse set 
of functions describes better the tails of the density and allows for a better (closer to the 
benchmark) coupling calculation. 

In contrast to the choice of NAKE, chosing the basis set has an effect on the accuracy 
of the calculated couplings. When opening the discussion to all the functional/basis set 
combination considered in this work, deviations are noticeable, especially for the QZ4P 
basis set for most of the systems at 3.50 A. For some dimers/XG functional combinations, 
there are deviations for all basis sets, see for example benzene dimer (BB) in Figure S3 
with the BLYP, PBEO, TPSS, and M06-HF. Thus, we distinguish two scenarios for these 
deviations, one is when one particular basis fails in conjunction with several XG functionals. 
We generally see this behavior for the more diffuse QZ4P set and almost never (a few outlier 
are the exception, such as the metaGGAs for AG and PBEO for BB) for the other sets. 
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Figure 4: MUE values as a function of the basis sets. PBE and PW91k are employed. See 
caption to Figure for additional details. 
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Figure 5: MRUE in % for the performance of each basis set. 


PBE and PW91k are employed. 
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The second situation is when two or all basis sets fail for a certain XC functional. This is 
attributed to a specihc shortcoming of the XC functional. 

In Figure SI the systems TH, PP, PY, IM and CD show an incorrect coupling at 3.50 
A (where the interaction between the two dyads is the highest) when the QZ4P set is used. 
Coupling values from 0.1 meV to 1000 meV are reported for these systems (see Table S2 
in the supplementary information) while the benchmarks are around 400 meV. The reason 
for this deviation is that FDE does not impose any constraint to a subsystem calculation to 
yield a charge localized electronic structure. The four factors responsible for the ability of 
FDE to yield charge-localized states are found to be systematic in their success of localizing 
the electronic structure on the subsystems. However, if the monomer basis set is large, 
one of the four factors becomes less effective, in turn increasing the chances of failure in 
the localization,^ or increasing the intersubsystem density overlap, in turn increasing the 
erroneous behavior of the NAKE.^ 


4.3 Effect of the XC functionals 


This section is devoted to the analysis of the performance of the XC functionals on the 
calculation of the electronic couplings with the FDE-ET method. Until now the basis set and 
NAKE functional correlations showed (besides the reported outliers) a relatively insensitive 
MRUE and MUE distribution along the considered range of distances. As we will see, 
this is not the case when varying the XC functional. The performance of each functional 
in all systems is presented in Figures ( |I0| ), ( pT| , and ( [I^ , where MUE (MRUE) is 

shown. From the hgures it is clear that all functionals behave well at the various distances 
for the majority of the systems. A different picture is presented for BB, in Figures and 
10 and for AC and BB, AA and TT in Figure S3. In these cases, even though the large 
majority of the XC/NAKE/basis set combinations are in good to excellent agreement with 
the benchmark results, some functional/basis set combination did not compare quantitatively 
to the benchmark. Looking at Figures and [T^ it is obvious that BB is the more problematic 
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Figure 6: MUE values as a function of GGA XG functionals. TZP and PW91k are employed. 
See caption to Figure]^ for additional details. 
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Figure 7: MUE values as a function of hybrid XC functionals. See caption of Figure for 
additional details. 
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Figure 8: MUE values as a function of the metaGGA and metahybrid functionals. See 
caption of Figure]^ for additional details. 
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Figure 9: Progressive polarization of the spin densities calculated for benzene (BB), naphtha¬ 
lene (NN), anthracene (AA), tetracene (TT) and pentacene (PP). Each column corresponds 
to a different XC functional. Benzene’s spin density is the most sensitive. 
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of the systems for BLYP. All other XC functionals perform well for this dimer. BLYP is 
surprisingly inaccurate at long and short ranges, with either small or large basis set. Although 
for this functional only, many possibilities were tried in order to improve the performance, 
for instance, additional freeze-and-thaw cycles were run (3 cycles are the default), a hner 
integration grid and more accurate density htting were also employed with no achieved 
improvements. This can be related to the known difficulties of semilocal functionals to 
model open shell systemsFor benzene and its derivatives, we notice that different 
functionals have an effect on the spin density polarization of the radical cation susbsystem, 
as exemplihed by the spin-density plots in Figure The BLYP functional produces the least 
spin polarized systems. A similar effect was reported for DNA nucleobase dimers^ where 
the spin density polarization was more pronounced for MP2 than for GGA XG functionals. 
In addition, we found the SGF to be slowly converging for BB, especially when BLYP is 
employed, as in the BB radical cation there is a degeneracy that is difficult to lift. Our claim 
is that this singular behavior of BLYP is related to its inability to match higher level of 
theory models for the radical cation subsystem, in turn undermining its ability to produce 
quantitatively correct couplings for BB. 

We do not provide here an explanation for the failure of metahybrids and metaGGAs for 
the AG system at long ranges. Use of these functionals in conjunction with FDE is so far 
untested and more investigations are needed to shed light on their behavior. 

In conclusion, GGA functionals are generally a good choice, with the unique exemption 
of BLYP for the benzene dimer. Table collects the best methods (i.e., the more stable 
for different system sizes) and PBE is reported as being the most accurate and transferable 
functional. The results for GGAs are in good agreement with the benchmark values, and 
in some cases they showed to be superior to non-local and MetaGGA functionals. Besides 
GGAs, B3LYP stands out as another valuable choice. 
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Figure 10: MRUE in % for the performance of GGA XG functionals. TZP and PW91k are 
employed. 


Table 2: Mean stistical values for the best XG-functional choices. 


Set 

MUE(meV) 

MRUE(%) 

MAX(meV) 

PBE/PW91k/TZP 

15.3 

7.1 

49.6 

PW91/PW91k/TZP 

18.0 

8.3 

48.4 

B3LYP/PW91k/TZP 

18.5 

8.0 

27.6 

M06-2X/PW91k/TZP 

29.1 

14.1 

90.0 
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Figure 11: MRUE in % for the performance of hybrid XC functionals. TZP and PW91k are 
employed. 


25 



















































S20 

10 


60 

50 

glO 

U30 

X 

S20 

10 


S20 

10 




Figure 12: MRUE in % for the performance of metaGGA and metahybrid XG functionals. 
TZP and PW91k are employed. 
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4.4 Electronic coupling dependence at different rotational angles 

In this section we focus on the sensitivity of the FDE method when the dimers are placed 
at different geometry conhgurations. The calculations were performed on ethylene and thio¬ 
phene dimers whose geometries were borrowed from Ref.^ The ethylene dimer was tested by 
rotating one ethylene molecule around the center of mass at 5 A intersubsystem separation. 
On the other hand, thiophene was classihed in three different type of rotations. First in 
a sandwich conhguration of the dimer, one thiophene molecule rotates around the center 
of mass, alternatively, the two thiophene molecules were rotated around the rotation axes 
that passes through the S atom, hnally, one thiophene rotates randomly around the center 
of mass. Intermolecular distances of 5 A, 6.75 A and 4 A were considered in order to keep 
a minimal distance between the hydrogen atoms of the dimer (for more details about these 
rotations on ethylene and thiophene see Figure 4 in Ref. IB7|1 . 

In Figure [T^ we report the performance of FDE-ET varying the XC functionals, as this 
was the one category that featured the largest deviations for the HABll previously studied. 
In the supplementary information, the analysis w.r.t. the basis set and NAKE functionals is 
also reported. Regarding the EE system (see Figure [T^), all functionals show appreciable 
deviations at 90° rotation angle. This is because the two double bonds are perpendicularly 
placed to each other and a nodal structure arises such that the overlap between the diabats 
is small. Because of this, numerical inaccuracies creep in the inversion of the transition 
overlap matrix to compute the transition density. Such a problem was detected already by 
us in a past study® and for which we proposed a solution based on the Penrose inversion 
of the transition overlap matrix. In this study, however, we purposely did not report values 
obtained adjusting this threshold. However, upon adjusting the threshold to a lower value, 
also the couplings at 90° compared quantitatively with the benchmark. In related methods, 
alternatives to the Penrose inversion have been proposed.l^^®^^®^ 

Generally, all functionals perform satisfactorily. The PBE and PW91 GGA functionals 
are in overall good agreement. However, there is a dependence on the basis set: the larger 
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Figure 13: Behavior of electronic coupling in thiophene dimer at (a) sandwich conhguration 
rotations, (b) simultaneous disrotation of the dimer and (c) random rotations, and (d) ethy¬ 
lene dimer. In each case the set of PBE functional (black circle), B3LYP (blue triangle) and 
M06-2X (magenta star) were used in conjunction with PW91k functional and TZP basis set. 
The orange line in (d) correspond to MRCI-I-Q. 

the basis set, the more accurate the coupling. In Figure S4, the basis set dependence of the 
electronic coupling with respect to the rotation angle is shown, although TZP and TZ2P 
perform well, the QZ4P basis set seems to yield the best results. Possibly because the 
distance between the ethylene monomers is larger than 4.5 A and thus in the long-range 
where we saw before the QZ4P performs best. 

The second test case is comprised of three different kinds of rotations of the thiophene 
dimer, see Figures [T^-[T^. In these three cases, the couplings are strongly dependent on 
the angle. The dependence is clearly due to the orbital overlap of the subsystem HOMOs, 
whenever there is a nodal structure (cancellation of different phases of the orbitals) the 
coupling will become negligible. By comparison of Figure 5 and Table IX on Ref. [671 with 
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the present results, the FDE-ET method ranks at the same level of CDFT with coupling 


values a bit lower than the reference. Figure l^i, shows that for the case considered, the 
coupling seem to be relatively independent form the chosen functional, or basis set (see 
Figure S4), and all functionals yield couplings within a few meV from the benchmark values. 


5 Conclusions 

The most important Ending in this work resides in the fact that GGA functionals coupled 
with a medium sized basis set and the PW91k NAKE functional allows the FDE-ET method 
to yield reliable electronic couplings as tested against high-level correlated wavefunction 
methods applied to an array of donor-acceptor dyads. We find the PBE functional to be 
the most transferable functional in each case considered having a MAX error lower than 30 
meV and an overall MRUE of a little over 7 %. This constitutes a success for the FDE-ET 
method. 

We analyze the performance of 10 XG functionals, ranging from GGAs to the Minnesota 
meta GGA functionals, and also hybrid functionals with Hartree-Fock exchange ranging 
from 10-30%, and metahybrid functionals with HF exchange in the 50-100% range. We 
extract from the statistics that the XG functionals are determinant in the performance of the 
electronic coupling. Gonversely, the NAKE functionals statistically do not play an important 
role (e.g., the couplings are relatively insensitive to their choice). In addition, our analysis 
of the basis set dependence shows that the QZ4P basis set (the largest set considered) is 
the most problematic as it often undermines the FDE convergence at short intermonomer 
separations - a problem already well documented in the FDE literature 

Overall, we show that by varying the three parameters considered in this study: XG 
functionals, NAKE functionals and basis sets, diabatic states are correctly generated with 
FDE-ET. In addition to the quality of the diabats, we provide convincing computational 
evidence that the FDE-ET method produces couplings which satisfactorily correlate with 
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the benchmark data. In conclusion, the FDE-ET is found to be a powerful tool for modeling 
CT (specihcally hole transfer) reactions. 
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